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Abstract 

nn , Epidemic models currently play a central role in our attempts to understand and control 

Q^ ' infectious diseases. Here, we derive a model for the diffusion limit of stochastic susceptible- 

Q . infectious-removed (SIR) epidemic dynamics on a heterogeneous network. Using this, we 

consider analytically the early asymptotic exponential growth phase of such epidemics, show- 
ing how the higher order moments of the network degree distribution enter into the stochastic 
behaviour of the epidemic. We find that the first three moments of the network degree distri- 
bution are needed to specify the variance in disease prevalence fully, meaning that the skewness 
s^ of the degree distribution affects the variance of the prevalence of infection. We compare these 

rf-\ . asymptotic results to simulation and find a close agreement for city-sized populations. 
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1 Introduction 



Epidemic models take a variety of mathematical form.s [2, 19], and are are now routinely used to 
inform policy on disease control and contribute towards public health plans [15, 34, 38, 4]. One 
of the most commonly investigated forms of models are susceptible- infectious-removed (SIR), in 
which everyone is either susceptible to the disease, infectious with it, or removed from the future 
C^ ' disease dynamics. This model is a simplification of reality, but provides a useful starting point for 

modelling diseases where previous infection confers long-lasting immunity, for example outbreaks 
of childhood diseases like measles, some respiratory illnesses like pandemic influenza, and historical 
pathogens such as smallpox [19]. 

Simple models of epidemics have assumed that all members of a population interact at a ho- 
mogeneous rate and therefore that at a given time every susceptible has an equal probability of 
contracting the disease from any infectious individual. A more realistic way to think of the spread 
of diseases can be that they spread through contacts between people, with these contacts describing 
a network of interactions. In the case of homogeneous mixing this network is under some conditions 
an Erdos-Renyi (ER) random graph. There have been many examples of using networks to study 
the spread of disease and two review papers [6, 10] compare several different approaches to network 
modelling. 



Longstanding generalisations of homogeneous mixing epidemic dynamics include heterogeneous 
network models [14] and regular clustered networks [18]. Regularity (each individual having exactly 
n links) is appropriate for some populations (and was originally motivated by spatially embedded 
ecological and veterinary applications) but for human respiratory transmission [29] and sexually 
transmitted infections [37, 17], this is almost definitely inaccurate and the contact network is highly 
degree-heterogeneous. In this work, we represent heterogeneity in degree using what has become 
a de facto standard: the configuration model [28]. 

As relevant theory has developed, through moment-closure approximations and other methods [33, 
6, 10, 23], the use of heterogeneous networks has become more and more common in the literature. 
Several areas of interest for epidemics have been studied. The invasion threshold of the infection has 
been calculated [12], and involves the mean and variance of the degree distribution. The final size of 
an epidemic with a given degree distribution along with the mean degree of the individuals infected 
has been derived [30]. In particular we note the use of a probability generating function [39, 26] 
to derive a small number of nonlinear ODEs that describe the dynamics of a SIR infection on a 
random heterogeneous network. In most cases where networks are used they are static, i.e. decided 
at one point in time and fixed that way rather than changing over time, though there are several 
examples where this assumption is not made [17, 40, 27]. There are also examples of networks on 
which the individuals have heterogeneous degree and are always changing contacts [32, 25], although 
there is some precedence for this modelling framework before the widespread use of networks as a 
conceptual tool [24] . Though the extremes of static or extremely dynamic networks are obviously 
unrealistic, they are useful to enable some analytic traction and also increase the ease with which 
we can simulate. More recently though there have been attempts to include some additional way 
of passing on the disease rather than just through the defined links of the network in an attempt 
to describe the fact that diseases can be passed between members of the population that do not 
have more than one interaction, i.e. a 'global' infection term [20]. 

Observed epidemics are noisy and unpredictable, which motivates the use of stochastic epidemic 
models [3]. While some asymptotic results for invasion and final size of stochastic epidemics on 
networks can be derived in a discrete-time branching process framework, if one is interested in 
transient dynamics then the natural model is an appropriate continuous time Markov chain. For 
an SIR-type epidemic model on an arbitrary graph with N nodes, such a Markov chain would 
involve 3^ ODEs, which quickly becomes computationally intractable. The method proposed 
by Ball and Neal [5] involves creating a configuration model network at the same time as the 
epidemic tree, and this stochastic process manifestly asymptotically converges on 2M ODEs in 
the deterministic large N limit, where M is the maximum number of contacts that the most well 
connected individual on a network has. Unfortunately, this can still be very large depending on 
the exact degree distribution and any attempt to cut down the number of equations by ignoring 
individuals who have more than a given number of neighbours, will inevitably ignore some of the 
most important individuals for the dynamics. Fortunately, recent work [11] has shown that a more 
sophisticated convergence proof leads to the smaller equation set of [39] . 

The desire to model stochasticity without a massive increase in dimensionality has led researchers 
to consider the diffusion limit. This general approach to stochastic processes is typically either 
attributed to van Kampen [31] or Kurtz [21, 22], but the basic idea is that provided that our 
population is sufficiently large, we can approximate the Markov process that describes the epidemic 
by a deterministic model together with appropriately scaled white noise processes that are defined 
by the transition rates between the states of the Markov chain. Such methods have been used to 



derive a low dimensional model in which properties of the noise in a stochastic epidemic model can 
be investigated analytically [1, 7], by Ross [36] to obtain expressions for the mean and variance 
of a metapopulation model, and by Colizza et al. [8] to model the effect of air travel on the 
spread of epidemics in a large-scale network. These models are attractive since they have the 
same dimensionality as the deterministic limit, but are actually stochastic, and corrections to the 
approximation are 0{N~^). 

In this paper we will apply the results of Kurtz to SIR-type epidemic dynamics on a configuration 
model network, as was done for SIS dynamics on a regular graph by Dangerfield et al. [9]. Using 
this we obtain a four-dimensional set of stochastic ODEs, from which we derive an analytical 
expression for the variance of the asymptotic early growth of an epidemic on a network given 
its degree distribution. We provide an argument that our approach is asymptotically exact, and 
simulate epidemics on various networks to confirm the utility of our analytical results. 



2 Network model 

The use of networks as a generalisation from homogeneous mixing is becoming one of the most 
widely used in epidemiological modelling. Contact between two individuals of the population we 
are considering forms a link between them. Once a link is established, the infection can be passed 
along it in either direction. Specifically what contact is represented depends on the disease, i.e. it is 
different for a respiratory infection compared to a sexually transmitted infection. These differences 
will affect how the contact network is constructed, but we use a common modelling framework 
once the network is known. 

The network can be described by a symmetric N x N matrix A , which has binary entries, Aij G 
{0, 1}, where the entry Aij will be 1 if nodes i and j have a link between them and otherwise. We 
do not allow self links in the network or multiple links between nodes. We make the simplifying 
assumptions that all links are of equal strength and that once a link is made it will be there 
throughout the spread of the epidemic. 

The degree distribution of a network, P(fc), specifies what the probability is that a node selected 
uniformly at random will have k neighbours. To construct an uncorrelated network with a given 
degree distribution, P{k), the configuration model is used [28]. The method for this is as fol- 
lows: 

• Each member of the population i is given a number of "half-links" or "stubs", ki, which 
is drawn from the degree distribution P{k). Once this is done we can define d^ to be the 
proportion of nodes in the network that have k neighbours. 

• Pairs of stubs are then picked uniformly at random and are then joined up forming an 
undirected edge between the two corresponding vertices. 

Later in the paper, when we simulate an epidemic on a network, this is the method that is used to 
create a network with the degree distribution that we require. Asymptotically, the differences due 
to different ways of dealing with repeated- and self-edges will be 0{N~^). Our choice for how to 
deal with these in a finite system is to obtain a list of all the nodes which have multiple connections 
between them and self-edges. One by one, the extra links between nodes will be broken, say between 
node i and j, and then a randomly selected and connected pair of nodes will also be selected, say 



nodes i' and j ■ We then connect i and i' together and j a-nd j together, which leaves the degree 
distribution unaltered. This is then repeated until there are no repeated- and self-edges. 



3 Diffusion model 

3.1 Model definition 

We assume a population of N individuals on a configuration model network. Individuals are 
stratified by their disease state S, /, or i?, and their degree on the network fc. Individuals of 
type Sfe become /fc at a rate equal to the product of the transmission rate r and their number of 
infectious neighbours. Individuals of type Ik become Rk at a rate 7. We are interested in the large 
N regime, in which we write \Sk\ for the expected number of susceptibles of degree fc, [7^] for the 
expected number of infectious individuals of degree /c, and \AB\ for the number of connected pairs 
of individuals on the network where one is type A and the other is type B. Omission of a subscript 
denotes implicit summation, e.g. \S\ = X^/ct'^fc]- 

To derive pairwise equations, we require knowledge of the neighbourhood of each node. This is 
due to the fact that when an infection or recovery takes place the number of pairs of type \AB\ 
will be changed in a way which is dependent on the neighbours that the node has. This is why 
it is not straightforward to write down a low-dimensional form for this process as the population 
size becomes large, since this requires the distribution of neighbours of each node. We therefore 
make the following assumption, that we will justify below: The distribution of neighbourhoods 
with X susceptibles and y infectious individuals around a susceptible individual of degree A;, is a 
multinomial with probabilities that do not depend on fc, i.e. 

<y^^ = {^^ (1 - P5 - qsf-''->lql , (1) 

\^i y / 

where, 

[SS] [SI] .^. 

Here the term {^ ) = k\/x\y\, is the multinomial coefficient. We have not detailed an assumption 
for the neighbourhoods of infected nodes here, as this requires more care to be taken. This will 
be detailed later, and will be denoted by D^ j,. With these assumptions, it is possible to reduce 
the state space of the Markov chain dramatically. In fact, a closed system is obtained that is of 
dimension four. The insight that allows this dimensional reduction is from [39], although our model 
uses pairwise notation and is derived explicitly from applying the results of Kurtz to an underlying 
stochastic process. We note that g{x) = J^k '^kX^ (where dk, as defined above, is the proportion of 
nodes with degree k) is the probability generating function (pgf) of the degree distribution. We also 
define and keep track of the variable 9 in our dynamical system. 9 is defined to be the proportion 
of degree one nodes that are still susceptible at time t. Infection down each link is assumed to be 
independent, which implies that the probability that a node of degree k being susceptible at time 
t will be 9^ . We can then write the number of degree k nodes that are still susceptible at time t 
as [Sk] ~ Ndk9{t)^ . If there are no degree one nodes in the population, then we can think of 9 
in terms of its relationship with degree k nodes. It is the 1/fc-th power of the probability that a 
randomly chosen degree k node is susceptible. 



We also note that many important features of the network can be written simply using the pgf. 
In particular: Ng{9) ~ N^j. dkO'^ = ^jJySk] = [S], and ^'(l) — X]fe ^"^^ ~ ^' where n is the mean 
number of links per node. 

To construct the deterministic limit of the epidemic on the large network (called the fluid limit by 
some), we must consider the change in the variables of our system due to an event. For the SIR 
model our events are infections of susceptibles and recoveries of infecteds. As we are considering 
the epidemic on a network, to fully describe the events, we are interested in how many neighbours 
someone has who are infectious and susceptible. Keeping this in mind we therefore consider two 
types of events. The first type of event is a susceptible of degree A;, who has x susceptible neighbours 
and y infected neighbours, becoming infected, i.e. going from S to /. The second is an infected 
individual of degree k, who has x susceptible neighbours, recovering from the infection, i.e. going 
from I to R. We denote these two events as er{k,x,y) and e-y{k,x) respectively, where r and 7 
are the rates of infection and recovery, and write £ for the set of such events. 

We then denote the rates of these events as /r(fc, x, y) and /^(fc, x). The rates are given by: 

/r(fc, X, y) = yT[Sk]Dly^^: , 
Mk,x)=j[h]Dl,. 

We then sum the changes in our variables due to an event, multiplied by the rate of that event, 
over all events to get the set of equations that we require. With this in mind, we have a state space 

p ^{e,[I],[SS],[SI]y obeying 

p - r ^ Ap, y [Sk] Diy^, + 7 E ^P^ i^k] Di^, , (4) 

k,x,y k,x 

where Ap^ and Ap^ are the changes to our variable under transmission and recovery respectively. 
We also note that as we are not keeping track of the variable [//] , we will not be concerned with 
the number of infected neighbours of an infective central node. Therefore when we are concerned 
with making an assumption about the neighbourhood of an infective, we will only be interested 
in the number of susceptible neighbours that it has. To work out the change in [SI] say due to 
an infection, we have to consider the neighbourhood of the infected node. The node in question 
was susceptible, becomes infected, meaning that the pairs which were [SS] pairs before infection, 
become [SI] pairs and the [SI] pairs will become [II] pairs. If our node had x susceptible neighbours 
and y infected neighbours, i.e. it formed x [SS] pairs and y [SI] pairs, then the change in [SI] will 
he X — y. We can do similar calculations for the other variables. We note that there is a subtlety 
in the calculations for as we require the node which becomes infected to be of degree 1, and it 
also depends on the number of nodes which are degree 1, given by Ndi. 

Doing this we obtain the following expressions for Api- and Ap^, 

Ap, = i-Sk^/Ndi, 1, -2x, x-yf , Ap^ = (0, -1, 0, -x)'^ , (5) 

which are the changes in the variables p from an infection and a recovery respectively. If we write 
down the exact but unclosed pairwise equation for our system p, then we get the following set of 
equations: 

Ng'{0) 
[I] = t[SI] - 7[/] (6) 

[SS] = -~2t[SSI] 
[SI] = t{[SSI] - [I SI] - [SI]) - -f[SI] . 



We note that the expressions that we obtain for the recovery terms (terms multiphed by 7) are 
the ones that we would require our assumption about the neighbours of infecteds to calculate from 
(4). Due to the fact that the equations concerning the recovery terms in (6) are exact and closed, 
we have two conditions that we require Z?^ j., to satisfy, viz. 



[h]Di^k^[I], 5]a;[4.]i?^,fc = [5/]. (7) 



/ , L f^ J X , /i, L J ' / , 



Now when we evaluate (4) using our main assumption (1) we get the model of [39] in pairwise 
notation. We note that the assumption of multinomially distributed neighbours allows us to get a 
relatively simple set of equations, as the summation over all events is equivalent to taking various 
moments of the distribution. We get 

[I]^t[SI]-j[I] 

g"{0) (8) 



[SS] = -2t[SS][SI]^^^,,^^^ 



9"{0) 



This can be thought of as a deterministic approximation to the underlying stochastic evolution of 
the epidemic and is the exact limit of the stochastic process, which we describe in the next section, 
when N -^ 00 and (1) holds. It is important to note that one cannot simply start with (8) and 
'add noise'. Defining appropriate events and rates, together with the explicit statement of (1), is 
essential. 

3.2 Diffusion limit 

The work of Kurtz [21, 22] also tells us that the noise around the deterministic limit given above 
should be a Gaussian centred around this limit with an appropriate density. These results require 
certain technical conditions to be satisfied in order to apply. In particular, they require a family 
of right-continuous, temporally homogeneous jump Markov processes with elements X^r indexed 
by an integer N. In our case, N is the number of nodes and the state of the Markov chain is 
i-^x y k'-^x k) where X^ f, is the number of degree-fc susceptible nodes with x susceptible and 
y infectious neighbours, and X^^. ^ is the number of degree-/c infectious nodes with x susceptible 
neighbours. We seek a limit 

TV ^00, Xly^^~^[Sk]Dly^^, Xlk~^[h]Dlk, (9) 

and then lump the equations derived to reduce the system dimensionality. The infinitesimal pa- 
rameters for the rate of going from state x to state x' also need to obey 

Ql^:^NF{N-^^,i^) , (10) 

for some function F . If we can write the rates in this form, the general results give us the exact 
formula to work out what these noise terms should be and the full stochastic equations in the 
diffusion limit are 



k'X.y k.x.y 

(11) 



-t- 7^ Ap^ [h] Dl^ + ^ Ap^ J7 [h] Di, ^7(0 



where Crit) and £,^{1) are independent standard Gaussian noise processes associated with trans- 
mission and recovery respectively. We note that there is no simple expression for the amplitudes 
of the noise processes, as the square root prevents us from taking the moments of the multino- 
mial distribution and then expressing these in terms of the variables of our system, as we do for 
the non-stochastic terms. Therefore we do not explicitly write the full stochastic system, as the 
equations would be in a similar (but more complicated looking) form to (11). 

We use methodology that is derived from [21, 22] and set out conveniently in [9], to analyse the 
variance in the infection levels during the early growth of the epidemic. We define the early growth 
period to be the time before there has been depletion to the pool of susceptibles which is significant 
enough to affect the rate of growth in the number of infecteds. Our starting point is assuming that 
the early growth of the infection is exponential, with rate r. That is, 

[/] = /exp(rf) , (12) 

where / is a constant related to the prevalence of the infection as the early asymptomatic behaviour 
commences. Note that an additional assumption in taking the diffusion limit is that this quantity 
should be significantly larger than 1 but also significantly smaller than the total population size 
N. We then use this to work out the early behaviour of the other variables, 

9=l-Kg[I]/N , 

[SS] - Ng'il) - Kss [I] , (13) 

[SI] - Ksi [I] . 

We now define the local covariance matrix associated with an event i.e. when a susceptible node 
becomes infected, or an infected node recovers. This is calculated as follows. 



^ij = Yl /e^P»,eApj,e , 



(14) 



ee£ 



where e £ £ means that e is either an infection or recovery event associated with a particular 
neighbourhood. Using the early growth Ansatz, specifically ignoring terms 0([/]^), we find that 
this can be written as G[J](t), where G is constant. We can now use the following equation from [9], 
which is again derived from the theoretical work of Kurtz [21, 22], 

ra"^ ~ B(T^ - cr^B^ = [Gexp(ri) - exp(Bi)Gexp(Bi)'^]/ , (15) 

where er^ is the time dependent covariance matrix of our state variables p, r is the early exponential 
growth rate, B is the Jacobian of the fluid limit of our system (8), evaluated using the early growth 
Ansatz (12) and (13), and G is given above. The aim is to find an expression for cr^, as this will 
give us the variance of the epidemic during the early growth phase. 

We now give specific details of how the various matrices in the above equation can be calculated. 
To calculate G, we calculate the matrix G and then divide by [/](i). To do this we substitute in 
the changes in variables Ap, which are given above by (5). This means that we can write G as 
follows, 

G = rJ2 y[Sk]Dly,kFr + 7 5^[4]/?^,feF^ , (16) 

k,x,y k,x 

where 
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(17) 



are the outer products of the change in each variable vectors Ap,- and Ap^ . 

Let us consider how this calculation is done in practice by working out the entry G^/j [551. We use 
the (2,3) entry of F,- = -2x and F-^ = 0, which gives us that G[/] [55] = -2t^^ ,j, ,^[S'fc] xy D^y,.. 
We can separate the sum over k from the sum over x and y, and we can use the fact that the (1,1) 
moment of a multinomial distribution with variables fc,p and q is given by fc(fc — l)pq. Therefore 
G[/]jSS] = — 2r ^j, ^ y[Sk] k{k — 1) psqs- We can then use our assumption about the probabilities 
of connecting to susceptibles or infecteds from (2) and sum over the indices. We get: 

G[jUss] = -2r[SS][SI]^^. (18) 

When we have calculated all the terms for the G matrix we will input (12) and (13) to obtain the 
correct matrix for the early growth period. We do this for evaluate this for G[/]_[ss] linearising 
it with respect to [I]/N, i.e. any term which is 0{{[I]/N)'^) will be ignored. We then obtain the 
following expression for G[/]^[55]: 



G[n,[5S] = -2T[l f '','^ + 0{[I]/Nf , (19) 



where we have used the fact that g'{9) and g"{9) will become ^'(1) and (?"(1) with the use of the 
early growth Ansatz, and then writing this in terms of g, where we define g*^"^ = g'"-'(l). This now 
gives us that, 

9"{g"-9' 
[g'Y 

We show all the entries of G in Appendix B. 



Gm,[SS, = -^r^-^^L-^ . (20) 



The B matrix from (15) as stated above is the Jacobian matrix of the deterministic limit of the 
system, evaluated using (12) and (13) and is therefore given by 



/ -^^jj, \ 
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r 



(21) 



Now we have calculated G and B we can solve (15) for cr^ 



3.3 Neighbourhoods around an infected node 

The discussion above has so far not dealt with the neighbourhood of an infective. This only 
becomes an important consideration for the bottom right term in F.^, where we have to work out 
the square of the change in [SI] pairs due to recovery of an infectious node; in particular we are 
interested in calculation of the quantity 

X:=^x2[4]Z?^_, , (22) 

and so the task is to determine the number of infectives of degree k, and the second moment of the 
distribution of the number of susceptibles around such infectives. Now, the rate at which infectives 
of degree k are produced is given by 

~ ^f*^"^ " "^'^'■^^' = -iVdfcA;00'=-i . (23) 



We also know that the probability that an infective of age a (a node which was infected length of 
time a ago) is still infective is given by e~^°' . We can therefore work out [Ik] at a given time t by 
evaluating the following integral: 

[h] = -Ndkk f e{t - a){6{t - a)f-^e-^'' da . (24) 

Now consider the neighbourhood around an such an infective. Ignoring the very first cases, every 
infectious individual must have been infected by someone, leaving fc — 1 individuals who are po- 
tentially susceptible. If the infection of the central node happened a time a ago, then each of the 
A; — 1 potentially susceptible neighbours has an independent probability e~'^°' of avoiding infection 
from that central node, and in the event that central infection is avoided and the neighbouring 
node is of degree I a probability of 9^^^ of avoiding infection from any other source. Summing over 
/ appropriately weighted then gives the general expression 

[h]Dlk - -Ndkk I 6{t - a){9{t - a))''-'e-''' Binfx fc - 1, ^^^jMe-A da , (25) 

where Bin(a;|fc,7r) is the binomial probability mass function, representing the probability that x 
of fc trials with independent probability of success tt will be successful. Since we are interested in 
early asymptotic behaviour, we then linearise this expression making use of the expressions 

6^1, e^-r^[I]. (26) 

This gives asymptotic results 

[4-]<fc^ [/]^ re-('^+^)''Bin(x|fc-l,e--'')da 

5 (1) Jo (27) 

,n^ + if 9" g'" ' 

=^x^[i]- ' 



g' \ T + 7' + 7 2r + 7' + 7 

We have now calculated the entirety of the G matrix and can express it in the form G[/], where 
G is constant, as is required to use the results of Kurtz. 



3.4 Rate of convergence 

We now consider the rate at which convergence to our model will happen. Our approach here is 
heuristic rather than formal, and is based on the fundamental contributions by Ball and Neal [5], 
and more recently Decreusefond et al. [11] We start by redefining the network and epidemic pro- 
cesses in ways that are less useful for practical calculation than our initial definitions, but which 
make the rate of convergence more clear. In each case, we start with N individuals indexed 
«,j, •■• 

In a configuration model process, we let individual i have Ki stubs, and start the network ajacency 
matrix Aij{t = 0) = 0,V7, j. Then the process is defined to take 

{Ki,Kj,Ai^j,A,^,) -> {K, - 1, Kj - l,Aij + 1, Aj,i + 1) at rate oc KiKj . (28) 

Running this process until the absorbing state Ki ~ 0,Vi gives the adjacency matrix of a con- 
figuration model network, although depending on what one decides about repeated- or self-edges, 
different corrections of 0{N^^) may arise [13]. 



In the epidemic process, individuals have non-independent random variables for their state Xi, that 
can take values S, I or R. Then 

{S^,Ij) -^ {h,Ij) at rate tA,j , (29) 

and infectious nodes recover at rate 7. The fundamental insight from Ball and Neal [5] is that 
these two processes can be combined. In this construction, every node in the network is given a 
number of half-links, which are then paired up as the epidemic progresses to form contacts between 
individuals in the following two ways: 

• An infected with I remaining half-links makes contacts at rate r/; if it links to a susceptible 
then the infection will be transmitted with probability 1. 

• When an infected recovers (which happens at rate 7) all of its remaining half-links will be 
independently paired off with remaining half-links in the population. 

The question of what a susceptible neighbourhood looks like at a given time in this picture is 
answered by halting the epidemic process, and running the configuration model process to com- 
pletion, then counting the neihbours of each susceptible. As noted by Decreusefond et Al. [11] 
(who use a slightly different construction but the same basic idea) at finite N this gives each sus- 
ceptible a multivariate hypergeometric neighbourhood, as they sample without replacement from 
the population. As N becomes large, this tends to a multinomial distribution with corrections 
O(iV-i). 

A similar argument can be made for the neighbourhood around an infectious individual, in terms of 
its links that remain unpaired. In the absence of susceptible depletion, it is then possible to make 
a stochastic version of the deterministic argument about the neighbourhood around an infective in 
§3.3 above. 

Since the diffusion limit deals with terms of 0(7V~^/^) and ignores terms of 0{N~-^), we are 
therefore justified in making the multinomial assumption for the neighourhood around a susceptible. 
We note that it therefore implies the exactness of several other deterministic approaches in the 
absence of clustering, such as the pairwise closure techniques from [18] and the maximum entropy 
moment closure method of Rogers [35] . 

3.5 Early growth variance 

As described above we can solve (15) for cr^. We did this through computer algebra, due to the 
complexity of expressions involved. The variance of the number of infecteds during early growth is 
shown in full in Appendix A. This can be simplified in the certain regimes. If we define the time at 
which the epidemic grows at the rate predicted in [12] as icariy, and the time at which the depletion 
of susceptibles affects the growth rate and we leave the early growth phase at tdcpiotod, then when 
the current time satisfies icariy <C i <IC idopictcd, we can simplify the expression in Appendix A. We 
note that this regime will exist in an infinite size network if the initial amount of infection in the 
network / is sufficiently small. It is dominated by a single term which involves the mean, variance 
and skew of the networks degree distribution along with the recovery and transmission rates of the 
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infection. In this limit, the mean and variance of prevalence obey 

Mean(/) -> /e""* , for r = r ( — - 1 ) - 7 



9' 
Ya.r{I) ^ .g'2(.g"' + (r + l)(.g" + .g')) 



(30) 



Mean(/)2 (^'2 - g"'^){{^ + T)g' - rg") 



where / is a constant related to the prevalence of infection as the early asymptotic behaviour 
commences at tcariy, and g'"-* = g^"''{\). 

Figure 1 shows how the variance in prevalence changes as skew in degree increases. Seeing that 
as the skew increases we get an increase in the variance of the epidemic during early growth, we 
would therefore see that if we had a network whose degree distribution is a power law, would show 
greater variation than a negative binomially distributed network with the same mean and variance 
would show. We think of this increase in variation being caused by the members of the population 
who are very well connected in the network. These people have been called super-spreaders in the 
past. If the disease reaches these people in the early growth phase of the epidemic then we can 
expect a rapid increase in disease prevalence, as there will be many [SI] pairs through which the 
disease can be passed, whilst if they do not get it in this early phase, there will not be this rapid 
increase, which generates the large variance in this stage of the epidemic. We note that for the / 
term in (30) we have no analytical traction on what this should be for a given epidemic. If we wish 
to compare this result to simulation, we will fit the value of / so that the analytical prediction and 
the simulated results agree at a given point. As this is simply multiplied by the other terms, fitting 
it simply scales the prediction up or down rather than fine tuning the prediction itself. 



4 Comparison with simulation 

We have also used simulation to test our conclusion that if we fix the mean and variance of the 
degree distribution, but have different skew values, then the network which has the higher value 
of skew will also exhibit a higher variance of infecteds during the early growth period. We wish to 
see whether we get a match between the analytical result that we have obtained and the simulated 
result. There are several difficulties to achieving this that we wish to note. First, the networks 
are extremely large in the case of the analytical results and we do not know how exactly fast the 
convergence to the answer should be (although the corrections are 0(7V~^) we do not know what 
prefactors multiply this term). We have run our simulations on networks of size 10^, representing 
a small city, which in the end seemed sufficient. 

Secondly, the analytical results only work for the early growth phase of the infection, which can 
be defined as the time when the pool of susceptibles is not significantly depleted. Whilst in an 
arbitrarily large network this can last for an arbitrarily long time, in a finite network, this is not 
true and will vary depending on the degree distribution of the network. This means that if we 
want to compare the early growth variance of two networks, we may have a very brief window in 
which we can actually do it. This is exacerbated by the fact that we will have to allow for a period 
at the very beginning of the epidemic to be discounted, as we want the system to approach some 
sort of early asymptotic behaviour unaffected by the random selection of initial infected. 

Finally, there were also several assumptions made for the analytical system such as how the number 
of triples in the network can be approximated by doubles and the pgf g{x). With the finite network 
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that we construct there is no guarantee that this wih be accurate; again, convergence is 0{N ■^) 
but we do not have knowledge of the prefactors. 

We implemented the Gillespie algorithm [16] to simulate the time evolution of the epidemic, on 
networks of size 10^. To obtain the correct mean early growth behaviour, we allow each simulation 
to achieve a certain number of infections (10^) and then we set the simulation time to zero and 
let the epidemic progress from there. By allowing this initial amount of infection we will be able 
to absorb the initial conditions that we choose and get to the average behaviour of the system. In 
the system of size 10^, allowing 10^ infections is also small enough that the susceptibles will not be 
depleted significantly enough to affect the rate of growth of the epidemic. We ran 10"^ simulations 
on two networks with the same mean and variance but with different skewness. We note that the 
same networks were used for each of the 10'^ simulations. 

Figure 3 shows the result of these simulations compared with the prediction given by (31) in 
Appendix A. We can see that as predicted the network with the higher skew exhibits more variance 
in number of infecteds in the early growth stage and also we see that the analytical prediction 
for what the variance should be, is a good fit for what the simulations show. The early time 
disagreement is as we would expect, since this is when the 0{N^^) corrections should be most 
significant. 



5 Discussion 

We have considered the spread of SIR- type infections on networks of heterogenous degree. Using 
the assumption that the neighbourhood around a node will have a distribution of susceptibles and 
infecteds which is multinomial (1), we have derived a low dimension deterministic approximation 
(8) to the stochastic dynamics of the infection which we can write down precisely in the form given 
by equation (11). A heuristic argument has been given to justify this assumption in the large N 
regime. 

Computer algebra was then used to calculate the covariance matrix of our model using an equation 
derived by Kurtz. This was then used to give us the variance of the early growth of the epidemic 
on the network and we extracted the dominant term in the ioariy ^ t ^ idcpictod regime which is 
given by (30). These analytic results were derived for networks of extremely large size. 

By comparing the solution that was derived analytically with simulations, we have been able to 
demonstrate that with a network of size 10^ we can get a strong agreement between the two as 
can be seen in Figure 3. 

We have shown that as the degree skewness of the network increases, so does the variance of the 
number of infecteds during the early growth. An implication of this result is that in the situation of 
an outbreak of a disease, if we are able to target the very well connected people in the network, then 
we can (as well as has been previously established) decrease the impact of the disease efficiently 
- but such an intervention would also mitigate the 'reasonable worst case scenario' in how many 
people will become infected by reducing the variability of the epidemic. 

Another potential application of our results is to improve the estimation of epidemiological parame- 
ters during early growth of an epidemic. The possibility of estimating the variability in prevalence, 
together with a relationship between such variability and underlying model parameters, could 
enhance statistical work on epidemic prevalence curves. 
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In conclusion, we have shown how a low-dimensional system of stochastic differential equations can 
be used to describe the diffusion limit of a stochastic epidemic on a heterogeneous network, and 
have drawn out epidemiological conclusions from this model. 
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Figure 1: Early asymptotic dependence of the standard deviation of infection prevalence, divided 
by infection prevalence, on the skew T of the network degree distribution. The curve is plotted 
from equation (30) in the main text. Parameter values are: mean degree « 5.4; variance in degree 
« 67.2; transmission rate r — 0.0308; recovery rate 7 = 0.1. Skewness of degree is varied between 
realistic values (0 and 100) to see how the variance of the asymptotic early prevalence of infection 
is affected. We see that as the skewness is increased we get a higher variability of prevalence. This 
is as expected, since the higher the skew, the more neighbours the most connected individuals of 
the population have, reducing the predictability of the epidemic due to chance events amongst this 
small but epidemiologically important group. 
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Figure 2: Test of the binomial assumption. 100 Monte Carlo realisations were run for parameters 
N = 10^, P(4) = 2/3, P{8) = 1/3, T = 0.8, 7=1. (a) shows a typical epidemic curve and four 
time points sampled, (b) shows the empirical histograms for / around an Sk node and 5" around 
an Ik node at each time point, and the equivalent binomial for susceptible central nodes, together 
with 95% prediction intervals across simulations (where each simulation is time-shifted to agree 
on the first time at which [/] — 100) to give an indication of finite-size effects. This shows the 
accuracy of the asymptotic results even at finite size. 
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Figure 3: Comparison of simulated results to analytical predictions. Dashed lines are simulations 

and full lines are analytical predictions. Simulations are on two different networks, which have 

the same mean and variance for their degree distribution (mean « 5.4 and variance « 67.2) but 

different skewness: 24.3 for red / black lines and 6.7 for the pink / grey lines. Transmission and 

recovery rates are t = 0.0308 and 7 = 0.1 respectively, (a) shows a period of time at which we 

have agreement in the growth of the number of infecteds between the two networks that is strongly 

in agreement with the theoretical prediction from [12]. (b) is also taken for this time and we can 

see that the theoretical prediction that we have described previously in this paper deviates slightly 

from simulation, with this most pronounced early on when 0{N~^) corrections should be most 

significant. 
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A Full form of Var(/) 

Here we reproduce the entire expression for the variance of the prevalence of the infection during 
the early growth phase. The variance of [/] at time t is denoted by a?jif^y 



[im g'2^(^j ^ j.){2j + r) 



7re-2T*(e*(^+'^) - 1 



+ (.9" - ,9') (7.9' + rig" g'))) + (7 + r)2Te^'{ 

g're-fHe'if+r) _ 1)(5"(^±^ + r) + @^ - g'r) 



gV(e*(^+'-) 


l)(a"( '^+'^ 1 t) I s"'(''+'') - 


-g'r) 


f 9'r{l- 


7 + r 


-g'r) 



'^ ^+-+^- -^ + e-^\g' - g")iW + rig" - 5')) 

+ (.9" - .9') (7.9' + rig" - .9'))) - 76-'"^* (75' + rig" - g')) (rig" - g')e'^^+^^ + 75' + g'r + g'r ~ g"n 

+ 5"7(7 + r)e''*(7 + r(^-l))) , 

(31) 

where g^""' = (?*-"^(l), 7 and r are the recovery and transmission rates respectively, r is the expo- 
nential rate of the mean early growth given in (30), and as in (30), / is a constant related to the 
prevalence of infection as the early asymptotic behaviour commences. 

B Full form of G 

We reproduce the matrix G here. Note that 7, r and g*^"^ are defined in the same way as in 
Appendix A and di is the proportion of nodes in the network which have degree 1. 

Gefi = r{g" - g')/NHW^ 
Gg,[^=G[iUe\^rig'-g")/Ng'\ 

Ge,[ss\ = G[ss],[e] = , 
Gg^^si]^G^siUe]^r{g" - g')/Ng'\ 

G[/].[/] = (7-r+^), 

G[iuss\ = G^ssui] = Mg' - g")g"/g'^ 

G[iusi] = G[sim = (.9' - 5") ((7 - r)g' + rg")/g'^ 

G[55],[55] = -4T(g'-5")(.9" + 5"')/5'^ 

G[ssusi] - 2r(.9' - g")g"'/g'^ 

G[si],isi] = ( - g'r + g"ir + (7 + r)/(7 + r + r))+ g"'{-f + r)/{-f + r + 2t)) / g' . 
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C Branching process approximation 

An alternative approaeh to modelling the early epidemic is to imbed a branching process in the full 
dynamics, which then allows us to generate results relating to the variance we would observe in the 
k-th. generation of infecteds. We let the degree distribution be D, so that the probability that a 
neighbour of yours has degree k is given by k¥{D = fc)/E[D], where P(-D = k) = dk and K{D) ~ n. 
Then we consider a Reed-Frost epidemic on the configuration model of the network given by our 
degree distribution. This is a discrete time model, meaning that we are concerned with the number 
of invectives in successive "generations" of the disease. We note that if you are an initial infective, 
then you are in generation of the epidemic, if you are infected by an initial infective, then you 
are in generation 1 and so on. For a Reed-Frost epidemic we denote the probability that infection 
will be spread across an edge from an infective to a susceptible by p. For example, if in generation 
n you consider an infected node, then there is an independent probability of p, that each of its 
neighbours will be an infective of in generation n + 1. 

Suppose that we have /„ infected individuals in generation n. Then in generation n + 1, there will 
be X^jll^ Bi infectives, where Bi is the number of 'children' of each infective i. Clearly the value 
Bi is dependent on the degree of infective i, which is say k. It is also clear that this will be at 
most k — I, as to be infective in the first place, one of the links to our node must have passed the 
infection down it, meaning that it has a non-susceptible neighbour at the other end of this link. 
In this construction, during the early growth period there will have been no significant susceptible 
depletion and therefore we will consider all the remaining fc — 1 nodes to still be susceptible. 

The Si's are independent and identically distributed. We want to calculate the probability gener- 
ating function for the B's. To do this we calculate, 

oo 

E(s^) = ^s''P(B = 6). (32) 

6=0 

To calculate F{B = b), we use the law of total probability as follows: 

oo 

F{B = 5) = ^ P(S = b\D = k)F{D = k), 

^='' (33) 

where we have the fc — 1 in the binomial coefficient and in the powers, to account for the fact that 
the fc-th link was responsible for the infection in the first place. 



So we have. 



ns'')^y^s'y(W\i-pt-'-'k '^' 



^" ^\ b r ^ '^' E(D)' 

fc=0 fe=o ^ ^ ^ ' 

1 v-^ v-^ /fc — 1 



^ ' b=0 k=0 



(34) 



EE \ M'(i-^)'""'^^^ 
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Now we swap the order of summation to get, 

CX3 / k—1 



^ ^ fe=0 

The expectation of B, is then given by taking the derivative of this expression with respect to s 
and then setting s to 1. When we do this we get, 

^/ s 1 „/ ^ E(D(D-l))p , , 

E B = ^-^g" 1 = , ' ■ 36) 

^ ' E{D)^ ^ ' E{D) ^ ' 

Similarly we get 

""'■^^) = EiD) + E(D) [ E(D) ) ■ ^^^) 

The expectation for the number of invectives in generation k follows 

E(4) = E ^ B, = E(B)E(/fe_i) = E(B)'=/o . (38) 

i=l 

The law of total variance then tells us that the variance of Ik is given by, 
var(/fe) = E(var( J^ B^\Ik-l)) + var(E J^ B^\h-l 



^ i=l 
Ik-l 



e[Y1 var(B,) j + var(E(Bi)/fe_i) , 



(39) 



= E(4_i)var(B) + var(/fe_i)E(B)2 , 
= /oE(B)'=-ivar(B) + var(/fe_i)E(B)2 . 

Using induction this gives us 

var(/fe) = IoE{B)''-'^{l +E{B) + . . . +E{B)''-^)vai{B) , 

,_i E(i3)'--l (40) 

-^°^^^) E(B)-1^'^^(^)- 

When we substitute in the the expression for the mean and variance of B from (36) and (37) 
respectively, we get, 

var(4)./or^^^^^^^V"ff^^^^^^V-l^ ' 



E(D) J W E(D) J J E(D(j-i))p -. 

E{D{D - 1){D - 2))p^ , E{D{D-l))p fE{D{D - l))pY\ ^^^^ 



E{D) E{D) V HD) 

■ ,, X A; — 1 / / ff \ ^ \ 1 / Iff 2 II / II \ '^ 

-I [l-E.\ (i^L-^] -1> I9P ,9P ( 9 P^ 



9' ) \\9' ) J ^-l\ 9' 9' \ 9' 

Now, letting k become large, and setting E{B)^ — e*"*, p — t/{t + 7), to correspond as closely as 
possible to continuous-time results, we obtain 

,r^ r 2rt9'9"'r^+9'9"r{T + j)-{g")^T^ 
g'{T + -/)[Tg" - (T + "f)g') 
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which is not in agreement with (30). This is as we would expect, because the reasoning in §3 above 
depends critically on the Markovian nature of the dynamics. In contrast, the branching process 
can account for non-Markovian dynamics, but does not allow for the fact that high-degree infective 
nodes create new infections more quickly than low-degree infective nodes. The two approaches are 
therefore best seen as complementary. 
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